function [z, df, d2f] = FUNC4(a,b,c,d,e,f)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This version is checked on Oct 30 2023
% By Jun Yu
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% a omega0_t
% b RL_t+1
% c H_t
% d RK_t+1
% e sigma0
% f sigma1

% checked
fun = @(x,y)exp(x+y).*exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)).*exp(-(y-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi));
ymax = @(x)log(b*c/(d*(1+c)))-x;

z = integral2(fun,a,100,-100,ymax);

% checked
funa =@(s)exp(a+s).*exp(-(a-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi)).*exp(-(s-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi));
da=-integral(funa,-100,log(b*c/(d*(1+c)))-a);

%checked
funb = @(x)c./(d.*(1+c)).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*exp(-(x-(-0.5.*e.^2)).^2./(2*e.^2))./(e.*sqrt(2.*pi));
db=integral(funb,a,100);

%checked
func = @(x)b./(d.*(1+c).^2).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi));
dc=integral(func,a,100);

%checked
fund = @(x)(c.*b)/(d.^2.*(1+c)).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*exp(-(x-(-0.5.*e.^2)).^2./(2*e.^2))./(e.*sqrt(2.*pi));
dd=-integral(fund,a,100);

% checked
fune = @(x,y)exp(x+y).*(-exp(-(x+0.5.*e.^2).^2./(2*e.^2))./(e.^2.*sqrt(2.*pi))+exp(-(x+0.5.*e.^2).^2./(2*e.^2))./(sqrt(2.*pi)).*(-((x+0.5.*e.^2).*e.^2-(x+0.5.*e.^2).^2)./e.^4)).*exp(-(y+0.5.*f.^2).^2/(2.*f.^2))./(f.*sqrt(2.*pi));
de = integral2(fune,a,100,-100,ymax);

% checked
fung = @(x,y)exp(x+y).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))./(e.*sqrt(2.*pi)).*(-exp(-(y+0.5.*f.^2).^2./(2*f.^2))./(f.^2.*sqrt(2.*pi))+exp(-(y+0.5.*f.^2).^2./(2*f.^2))./(sqrt(2.*pi)).*(-((y+0.5.*f.^2).*f.^2-(y+0.5.*f.^2).^2)./f.^4));
dg = integral2(fung,a,100,-100,ymax);

df=[da db dc dd de dg];


% second order derivatives
% a omega0_t
% b RL_t+1
% c H_t
% d RK_t+1
% e sigma0
% f sigma1
    % define the probablity density function
    function [f0] = desity(x,sigma_under)
      f0 = 1/(sigma_under*sqrt(2*pi))*exp(-(x+0.5*sigma_under^2)^2/(2*sigma_under^2));
    end
    
    % call the probability density function
    f0_omega0bar = desity(a,e);
    f1_omega0bar = desity(log(b*c/(d*(1+c)))-a,f);
    
    % first order derivative of probability density function with respect
    % to sigma
    function [f0prime] = desityprime(x,sigma_under)
      f0prime = -1/(sigma_under.^2*sqrt(2*pi))*exp(-(x+0.5*sigma_under.^2).^2/(2*sigma_under.^2))+ 1/(sqrt(2*pi))*exp(-(x+0.5*sigma_under.^2).^2/(2*sigma_under.^2))*(-(((x+0.5*sigma_under.^2)*sigma_under.^2 -(x+0.5*sigma_under.^2).^2)/(sigma_under.^4)));
    end
    
    % Second order derivative of probability density function
    function [f0prime2] = desityprime2(x,sigma_under)
      f0prime2 = 2/(sigma_under.^3.*sqrt(2.*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)) ...
               + 1/(sqrt(2*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)).*((((x+0.5.*sigma_under.^2).*sigma_under.^2 -(x+0.5.*sigma_under.^2).^2)/(sigma_under.^5)))...
               + 1/(sqrt(2*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)).*(-((2.*x.*sigma_under.^2+2.*sigma_under.^4-6.*sigma_under.^2.*(x+0.5.*sigma_under.^2)+4.*(x+0.5.*sigma_under.^2).^2)./(sigma_under.^5)))...
               + 1/(sqrt(2*pi)).*exp(-(x+0.5.*sigma_under.^2).^2./(2.*sigma_under.^2)).*(-(((x+0.5.*sigma_under.^2).*sigma_under.^2 -(x+0.5.*sigma_under.^2).^2)./(sigma_under.^4))).*(-(((x+0.5*sigma_under.^2).*sigma_under.^2 -(x+0.5*sigma_under.^2).^2)/(sigma_under.^3)));
    end
    f0prime_omega0bar = exp(-(a-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))*(-(a-(-0.5.*e.^2))./(e.^2));

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For da functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    funaa = @(x)exp(a+x).*(f0_omega0bar+f0prime_omega0bar).*exp(-(x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi));
    daa1=integral(funaa,-100,log(b*c/(d*(1+c)))-a);
    
    % get daa %checked
    daa = (b*c/(d*(1+c))).*f0_omega0bar.*f1_omega0bar-daa1;
    
    % get dab %checked
    dab = -(b*c/(d*(1+c))).*f0_omega0bar.*f1_omega0bar./b;
    
    % get dac %checked
    dac = -(b*c/(d*(1+c))).*f0_omega0bar.*f1_omega0bar./(c.*(1+c));
    
    % get dad %checked
    dad = -(b*c/(d*(1+c))).*f0_omega0bar.*f1_omega0bar./(-d);
    
    % get dae %checked
    f0prime_omega0bar_sigma0 = desityprime(a,e);
    funae = @(x)exp(a+x).*f0prime_omega0bar_sigma0.*exp(-(x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi));
    dae = -integral(funae,-100,log(b*c/(d*(1+c)))-a);
    
    % get dag %checked
    funag = @(x)exp(a+x).*f0_omega0bar.*(-1/(f.^2.*sqrt(2.*pi)).*exp(-(x+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((x+0.5.*f.^2).*f.^2 -(x+0.5.*f.^2).^2)./(f.^4))));
    dag = -integral(funag,-100,log(b*c/(d*(1+c)))-a);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For db functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % get dba %checked
    dba = dab;
    
    % get dbb %checked
    funbb = @(x)(c/(d*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./b;
    dbb = integral(funbb,a,100);
    
    % get dbc %checked
    funbc1 = @(x)(c/(d*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(c.*(1+c));
    funbc2 = @(x)(1/(d*(1+c)^2)).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)));
    
    dbc = integral(funbc1,a,100)+integral(funbc2,a,100);
    
    % get dbd %checked
    funbd1 = @(x)(c/(d*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(-d);
    funbd2 = @(x)(c/(d^2*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)));
    
    dbd = integral(funbd1,a,100)-integral(funbd2,a,100);
       
    % get dbe %checked
    funbe = @(x)(c./(d.*(1+c))).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-1/(e.^2.*sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4))));
    dbe = integral(funbe,a,100);
    
    % get dbf %checked
    funbg = @(x)(c/(d.*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(-1/(f.^2.*sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2/(2.*f.^2))+ 1/(sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2/(2.*f.^2)).*(-(((log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).*f.^2 -(log(b.*c./(d.*(1+c)))-x+0.5.*f.^2).^2)./(f.^4))));
    dbg = integral(funbg,a,100);
      
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For dc functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
    % get dca %checked
    dca = dac;
    
    % get dcb %checked
    dcb = dbc;
    
    % get dcc; %checked, a problem found
    funcc1 = @(x)(b/(d*(1+c)^2)).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(c.*(1+c));
    funcc2 = @(x)(2*b/(d*(1+c)^3)).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)));

    dcc =  integral(funcc1,a,100) - integral(funcc2,a,100);
    
    % get dcd %checked
    funcd1 = @(x)(b/(d*(1+c)^2)).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2))./(f.^2)))./(-d);
    funcd2 = @(x)(b/(d^2*(1+c)^2)).*(exp(-(x-(-0.5.*e.^2)).^2/(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)));

    dcd = integral(funcd1,a,100) - integral(funcd2,a,100);
    
    % get dce %checked
    funce = @(x)(b/(d.*(1+c).^2)).*exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2/(2.*f.^2))./(f.*sqrt(2.*pi)).*(-1/(e.^2.*sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4))));
    dce = integral(funce,a,100);
    
    % get dcg %checked
    funcg = @(x)(b/(d.*(1+c).^2)).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(-1./(f.^2.*sqrt(2*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).*f.^2 -(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2)./(f.^4))));
    dcg = integral(funcg,a,100);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For dd functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    
    % get dda
    dda = dad; %checked
    
    % get ddb
    ddb = dbd; %checked
    
    % get ddc;
    ddc = dcd; %checked
    
    % get ddd %checked
    fundd1 = @(x)(c.*b./(d.^2.*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2))./(f.^2)))./(-d);
    fundd2 = @(x)(2.*c.*b/(d.^3.*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(exp(-(log(b.*c./(d.*(1.+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)));

    ddd = -integral(fundd1,a,100) + integral(fundd2,a,100); 
    
    % get dde %checked
    funde = @(x)(b.*c./(d.^2.*(1+c))).*exp(-(log(b.*c./(d.*(1+c)))-x-(-0.5.*f.^2)).^2./(2.*f.^2))./(f.*sqrt(2.*pi)).*(-1./(e.^2.*sqrt(2*pi)).*exp(-(x+0.5.*e.^2).^2/(2.*e.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)/(e.^4))));
    dde = -integral(funde,a,100);
    
    % get ddg %checked
    fundg = @(x)(b.*c./(d.^2.*(1+c))).*(exp(-(x-(-0.5.*e.^2)).^2./(2.*e.^2))./(e.*sqrt(2.*pi))).*(-1./(f.^2.*sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2./(2.*f.^2))+ 1./(sqrt(2.*pi)).*exp(-(log(b.*c./(d.*(1+c)))-x+0.5.*f.^2).^2./(2.*f.^2)).*(-(((log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).*f.^2 -(log(b.*c./(d.*(1.+c)))-x+0.5.*f.^2).^2)./(f.^4))));
    ddg = -integral(fundg,a,100);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For de functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    
    % get dea
    dea = dae; %checked
    
    % get deb
    deb = dbe; %checked
    
    % get dec
    dec = dce; %checked
    
    % get ded
    ded = dde; %checked
    
    % get dee; %checked
    % may be problematic
    funee = @(x,y)exp(x+y).*(desityprime2(x,e)).*exp(-(y+0.5.*f.^2).^2./(2.*f.^2))./(f.*sqrt(2.*pi));
    ymax = @(x)log(b.*c./(d.*(1+c)))-x;

    dee = integral2(funee,a,100,-100,ymax); 

    % get deg %checked
    funeg = @(x,y)exp(x+y).*(-1./(e.^2.*sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))+ 1./(sqrt(2.*pi)).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2)).*(-(((x+0.5.*e.^2).*e.^2 -(x+0.5.*e.^2).^2)./(e.^4)))).*(-1./(f.^2.*sqrt(2.*pi)).*exp(-(y+0.5.*f.^2).^2./(2.*f.^2))+ 1/(sqrt(2.*pi)).*exp(-(y+0.5.*f.^2).^2./(2.*f.^2)).*(-(((y+0.5.*f.^2).*f.^2 -(y+0.5.*f.^2).^2)./(f.^4))));
    deg = integral2(funeg,a,100,-100,ymax);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % For dg functions
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    
    % get dga
    dga = dag; %checked
    
    % get dgb
    dgb = dbg; %checked
    
    % get dgc
    dgc = dcg; %checked
    
    % get dgd
    dgd = ddg; %checked
    
    % get dge
    dge = deg; %checked
    
    % get dgg %checked
    % may be problematic 
    fungg = @(x,y)exp(x+y).*exp(-(x+0.5.*e.^2).^2./(2.*e.^2))./(e.*sqrt(2.*pi)).*(desityprime2(y,f));
    ymax = @(x)log(b.*c./(d.*(1+c)))-x;

    dgg = integral2(fungg,a,100,-100,ymax);   
    
    
d2f = [daa dab dac dad dae dag; 
       dba dbb dbc dbd dbe dbg;
       dca dcb dcc dcd dce dcg;
       dda ddb ddc ddd dde ddg;
       dea deb dec ded dee deg;
       dga dgb dgc dgd dge dgg];


end



